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We extend the Slepyan solution of the problem of a steady-state crack in an infinite ideally 
brittle lattice model to include dissipation in the form of Kelvin viscosity. As a demonstration of 
this technique, based on the Wiener-Hopf method, we apply the method to mode III cracks in a 
square lattice. We use this solution to find the critical velocity at which the steady-state solution 
becomes inconsistent due to additional bond-breaking; this point signaling the onset of complex 
dynamical behavior. 
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I. INTRODUCTION 



Recent experiments on dynamical fracture (Fineberg, et al., 1991, 1992| ) have shown that cracks in brittle materials 
become unstable above a critical velocity. This instability is associated with the formation of a roughened fracture 
s urface, with an increased d issipation of elastic energy, and with complicated velocity dynamics - for a review, see 
( Fineberg and Marder, 1999]). Alth ough there are hints of such an instability in the traditional continuum form ulations 
of ideally brittle crack s ( |Yoffc, 1951 ), a systematic treatment doe s not appear possib le. Indeed, recent attempts ( Langer 
and Lobkovsky, 1998| ) to utilize a cohesive stress modification (Barcnblatt, 1959) of the continuum elastic equations 
in order to address this problem have been rather unsuccessful. 

In the absence of a compelling continu um approach, one mu s t deal in some manne r with discrete dynamics on the 
"atomic scale". In this regard, Slepyan ( Kulamekhtova, 1984; [Slepyan, 1981 , 1982) pioneered the idea of studying 
lattice models in which atoms interact (via piece-wise linear springs) with their neighbors i n a predetermined lattice 



geometry. This conceptual framework was further developed by Marder and collaborators (Marder and Gross, 1995 
Marder and Liu, 1993| ). Here, one can directly obtain the relationship between the crack tip velocity and the imposed 



driving and furthermore one can search for conditions which do not allow stable steady crack motion. Clearly, 
lattice models are not fully realistic, especially at displacements that are not small compared to the underlying 
lattice spacin g; to be more realistic, one must resort to mole c ular dynamics simulations including all possible atomic 



1995 



inter a ctions (|Abraham, et al., 1994| |Gumbsch, ct al., 1997| ; |Holland and Marder, 1997| , |1998| ; Pmcltchcnko, et al 



1997; Zhou, ct al., 1997). However, the analytic tractability of these models as well as indications (Marder and Gross 



Pla, et al., 199S ; Bander and Ghasias, 199S ) that they do contain the essential mechanism responsible for the 



observed dynamical instability make them well worthy of serious attention. 

Slepyan recognized that in an ideally brittle material (where each spring is linear until a displacement at which 
it completely breaks) one can use Fourier methods to find the lattice analog of a traveling wave solution. In this 
solution, each point on the lattice (at the same transverse coordinate) undergoes the same time history of motion 
as any other, albeit with some time delay. Once one obtains the solution, one must check that all springs assumed 
to be linear have displacements that are in fact below the breaking threshold. The violation of this assumption 
by the steadily propagating solution s ignals the onset of more complex dynamical behavio r, at least in qualitative 
accord to what is seen experimentally (Fineberg and Marder, 1999; Marder and Gross, 1995). Using the Wiener-Hopf 
technique, Slepyan solved for the velocity-driving curve in the limit where the width of the lattice transverse to the 
crack direction goes to infinity. 
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One important consideration in all models of brittle fracture concerns dissipation mechanisms. From our perspective, 
it makes sense to put in dissipation at the lattice scale in such a way so as not dominate the large-scale continuum 
elastic field. If in addition one demands an equation that is local in time, one is led ( Langer, 1992j ; Pla, et al., 199S ) 
to the introduction of a Kelvin viscosity which dissipates energy proportional to the rate of change of spring lengths. 
In the naive continuum limit, this gives rise to a third derivative (two space, one time) term. That means that if the 
viscosity is chosen to be 0(1) on the lattice scale, i t will scale to zero as far as the macrosco pic dynamics is concerned. 
We saw this explicitly in a previous set of papers (Kessler, 2000; Kessler and Levine, 1998) which solved this problem 
for finite width lattices and considered the nature of the solution as the width became large. Nevertheless, the viscosity 
can have a considerable effect on aspects of the solution which depend on the microscopic details, namely the crack 
speed (for a fixed stress intensity factor) and the self-consistency of the traveling wave ansatz. 

Th e rest of this paper is organized as fo llows. In the next section, we formally define the model, employing Slepyan's 
idea ( Kulamckhtova, 1984; Blepyan, 1981) of replacing the driving at some external boundary with some local forcing 
on the crack surface. This allows for the solution in Section III of the strictly infinite lattice model; the connection 
to a finite problem with some displacement driving is made by appealing to the universality of the microscopic crack 
solution given a fixed stress intensity factor. In the following section, we discuss the numerical evaluation of the 
velocity curve from its formal expression. Following that, we turn to the aforementioned self-consistency condition 
and study the effect of dissipation on the critical velocity. We conclude with some observations in the final section. 



II. SQUARE LATTICE MODEL 

We wish to study mode III cracks and the effects thereupon of Kelvin viscosity. Our model thus consists of a 
square lattice of mass points undergoing (scalar) displacements out of the plane. The lattice extends infinitely long 
in both the x (along-crack) and y (transverse) directions. The lattice points are connected to their nearest neighbors 
by ideally brittle springs which behave linearly with spring constant 1 until some threshold elongation 2e at which 
point they irreversibly crack. All the (uncracked) springs have a viscous damping 77. Let us assume that the crack 
corresponds to a sequential breaking of the vertical bonds between the mass points at rows y = 1 and y = 0. The 
equation of motion for the masses in rows y > 1 reads 



l(x, y) = (l + r)^\ (u(x + l,y) + u(x -1,2/) + u(x, y + 1) + u(x, y - 1) - 4u(x, y)"j 



(1) 



We allow the vertical bonds along the crack surface to have a different spring constant k and damping parameter fj, to 
model the effect of having weak links along the crack surface. The equation for y = 1 then reads, using the assumed 
symmetry u(x, y) = —u(x, 1 — y), 

u(x,y) = (l +11^ (u(x + l,l) + u(x- 1,1) +tt(x,2) - 3u(x, lj] 



+8(e - u(x, 1)) (k + fj^j ( - 2u(x, 1)) 



(2) 



It is easily verified that the equations for y < are consistent with the above symmetry. We can, without loss of 
generality choose e = 1/2 (so that the threshold elongation is unity). Note that in these units, the elastic wave speed 
is unity, so all velocities are dimensionless, expressed as fractions of the wave speed. 

We are interested in steady-state cracks, described by the Slepyan traveling wave ansatz, 

u(x,y,t) = u(x-vt,y) (3) 

which implies that every mass point in a given row undergoes the same time history, translated in time. We choose 
the origin of time such that it(0, 1) = e, so that it represents the moment of cracking of the spring in column 0. If we 
define r = x — vt, the equations of motion for y > 1 become 

( 1_7?,; J^) (u(T + l,y) + u(T-l,y)+u(T,y-l) 

+ u(r, y-1)- 4u(r, y)) - v 2 ^^ - . (4) 
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For y = 1, we separate out the terms proportional to 6{—t), giving 

d 



1 — rjv 



^ (u(t + 1,1) +u(t- 1,1) + u(t, 2) - 3u(t, 1)) 



The driving term for this equation has the form 



F„ - 2 [k-rjv-?- 



u(r,l) 



(5) 



(6) 



In this expression, we have inserted an external force -Fb which acts on the crack boundaries. This force is an artifice 
which serves as an inhomogeneous source term, in place of introducing the driving through the boundary condition 
on the top and bottom surfaces. We will see later how to choose this function. 

We Fourier transform with respect to r (see the appendix for our conventions regarding Fourier transforms and 
Fourier integrals). If we assume that u{q,y) — Xg[£(9)P _1 f° r V ^ 1; we obtain from Eqn. (0) the dispersion 
relationship 



(1 + iqr 1 v){e- 1 ' 1 + + £ + 1 - 4) + q 2 v 2 = 



This can be rewritten as 



1/2 



£1/2 



) 2 =4sin 2 ^- 



2 2 

q v 
1 + zqryw 



(7) 



(8) 



which therefore implies that £ 1/2 = ±h/2 ± y / {h/2) 2 + 1. If we define 



- 2 = h 2 + 4 



then we can write 



_ (r - h) 2 _ i h(h-r) _r-h 



(9) 



Note that we have must choose the sign that corresponds to a decaying solution in the y direction, i.e. |£| < 1. 
If we substitute this ansatz into Eqn. ^ , we get 

— 2(fc + iqfjv) + (1 + iqrjv)(l — x = -& F 



where a F is the Fourier transform of a. Noting that —1 

Xq = 



1 _ 

t 

2a F 
N(q) 



r{r+h) 
2 ' 



we obtain 



with 



N(q) = h(r + h)(l + ir/vq) + A(k + ifjvq) 



(10) 



(11) 



(12) 



This equation is implicit, since the driving term a F on the right hand side still depends on the displacement field 
u(t, 1). In the next section, we solve this equation via the Wiener- Hopf technique. 



III. WIENER-HOPF SOLUTION 

Let us define V(t) = 2u(r, 1), the bond elongation between rows y — and y = 1. In the appendix, we define V + (q) 
and V - (<7) as the decomposition of V(r) into terms which are analytic in the upper and lower half planes respectively. 
This allows us to write Eq. (0) as 



a F (q) = F-(q) + (k + ifjvq) V~ (q) - ijv V (0) 



(13) 
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where we have employed a similar breakup for the driving term Fq, and where we have utilized Eq. (|62j). 

Substituting the formula for a F into Eqn. ( |To| ) for \ = (V + + V~)/2 leads to a closed form equation for V; namely 



Defining S by 



we can rewrite Eq. (|l4|) as 



, _ 4(fc + iijvq) AF^ _ Arjv V (0) 

' { N(q) ' N(q) N(q) 

h(l + irjvq) 



S=l- 



4(fc + irjvq) 



5 V" 



)]V 



h{\ + i7]vq) + (r — h)(k + ifjvq) ' 
(1 - S)F Q 



-(l-S)V(O) 



iqrjv 



iq-qv 



(14) 



(15) 



(16) 



We will solve this equation using the Wiener-Hopf technique. The key is to factor S into a product of two pieces, 
S(q) — S + (q)S~(q), each of which is regular in the upper (S + ) or lower (S~) half plane. Doing so, we get 



V+ nv ( 1 

7TX + S~ V" +- '— — — - S 1 " 



V(0) = 



1 



F- 



S+ k 



iqrjv 



We then have to deal with the extra singularity at q = We do this by a simple subtr- 



action 



S+ k + iqr/v \ S+ 



1 



l-S Fn 



5^ 



(17) 



(18) 



We have now to choose Fq, our external forcing. In taking the width to infinity, in essence we are solving for the 
"inner" solution in the sense of boundary-layer theory; i.e., only on the scale of the lattice spacing. The external 
forcing in the finite width problem only acts on the large scale, and does not vary on the lattice scale. Thus, F must 
have support only at q — 0. Following Slepyan, then, we choose Fq such that 



1 



S + (k + iqfjv) 



F- = 2tt.B%) = B 



1 



1 



iq + -iq + 



where B is a constant. This then gives us the the two separate equations 



= - 



fjv V (0) 
k + iqfjv 

rjv V (0) 



1 - 



S+ 



1 



H rjv 



BS+ 
-iq + 



B 



{iq + 0)S- 



k + iqrjv \ S~S + \ ik_ 

As shown in Eqs. (|5|) and ( ]6l| ) we can solve for V(0) from Eqs. ( |20a| ) or ( |20b) ) by multiplying the first by — iq 
or the second by iq — > +00. This yields 

V(0) = BS + |, ■ 



Substituting this back into Eqns. (20), and defining r/k = fj/k, we find 

, BS+ r] k vB r 

V + = S 

(~iq + 0)(l + iqr] k v) 1 + iqrj k v 



(19) 

(20a) 

(20b) 

-» +00 
(21) 

(22a) 



V- = t ^ ^ + ^^-^ + | * S- (22b) 

(iq + 0j(l + iqr]kV)b 1 + iqrjkv q ~ n k « 

These equations directly determine the displacements given the strength of the external driving, B. To get a 
physical handle on the meaning of this constant, we note that Eqns. (g0|) can be directly solved for the leading 
behavior of V ± as q goes to 0, using the fact that in this limit, 

S ± ~K kX VA(q±iQ) 1 ' 2 (23) 



5 



rr— 2" 

where A — 2 T , and k is a constant, reflecting the arbitrariness in how we perform the decomposition of S. We 
find 

v + - ^ ( _; g +o)i/ 2 ( 24a ) 

v -7I(^Toj^ (24b) 

This is just the expected (Fourier transform of) the singular solution of continuum elasticity as one approaches the 
crack tip. Thus, as in all boundary-layer problems, the large-distance limit of the "inner" solution matches on to 
the short-distance limit of the "outer" solution, confirming the justification of the choice of Fq above. We need to 
relate the strength of the short-distance singularity of the continuum solution to the driving displacement A. This 
can be done directly by solving the continuum elastic equations. It is easier, however, to follow Slepyan and derive 
this relation by energy considerations. One can calculate (Slepyan, 1981) the flux of energy into the boundary-layer, 
or "process" zone, obtaining 

ikn 2 B 2 , , 

T=^—v (25) 

Now consider the flux of energy To being used to break bonds, which is simply related to V(0) 

k(B S 



To = -± ^ VkV/ v . (26) 

Then, using the idea that the Griffith's displacement is determined by exactly the condition that all the energy flux 
is needed just to break the bonds, we have 

^- (27) 

1=^ 



A G V T S+ 



Notice that the factor of k guarantees that the result is invariant with respect to how the decomposition of S is 
performed. 

In order to make contact with other treatments, we consider the standard case k = 1 and rjk = fj — f], in which all 
bonds are equivalent. In this case, the function S has a particularly simple form, 

S( q ) = ~ r = | (28) 

with H 2 = 4(1 + iqryri) sin 2 | — q 2 v 2 and R 2 = 4(1 + iqrjv) (sin 2 | + 1) — q 2 v 2 . We can represent these functions in 
term of the products over all their roots 

1/2 / \ 1/2 

9 /, q 



I! .ill 



h = + ioy/ 2 d - *o) 1/2 n ( 1 - ^rj y 1 ^ ) (29) 



^ = 2 n( i -4) (i--) (30) 

where Qni^t) are the (nonzero) roots of H 2 (R 2 ) in the lower half plane and Qniln) are corresponding roots of 
H 2 (R 2 ) in the upper half plane. This decomposition allows us to write down explicitly the factors S + , S~: 



g+(g)= (1 "^ )1/4 (9+i0)1/2 n(T 



1 - 



_2_\ V2 

Qi \ 



Q 
„+ 



1 _ _q_\ 1/2 



(31a) 
(31b) 
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Examining the small q behavior, we find the k — 1, so that substitution of S + into Eqn. ( p7j ) gives 

i \ 1/2 



A 



(2r/v) 



1/2 / 1 



>/<'</„ 



If we explicitly take out of the product in Eq. ( |32|) the one imaginary q + , which we will denote as q^ , we get 

1/2 



(32) 



1 



A G (l-W 2 )l/4 



! 2(l + iv m +)_ n+iv m +Q; 

n^O.m 



ivr/Qr. 



(33) 



This formula was obtained by Kessler (2000), who started from a finite lattice of transverse size N y and considered 
the N y — > oo limit. Note that this can be done for k = 1 on a square lattice, (the case being considered h ere), but not 
for gener al k or for modes I or III on a triangular lattice; the method here works for these cases as well (Pechenik, et 
al., fcOOOP. 

Next, it is easy to see that in the limit r\ — > + , we have S\ q= j_ — 1 and that therefore 



A 

Ag 



s- 



1/2 



;i/2 



n 



QmQn 



(34) 



Furthermore, in this limit q$ 



-q Q . All other complex roots cancel out, leaving only real roots. Hence, 
Ac" 



n 



over real 



Qm Qn 



(35) 



a result originally obtained by Slepyan. Note that in our treatment, the presence of one of these real roots in either 
the numerator or the denominator of the product expression depends entirely on the half-plane from where it came 
as we took the limit 77 — i- + . To see what this condition implies, let us denote by H 2 (ot R 2 ) as /; then the small 
imaginary part of a specific real root, i5, satisfies the equation 



iS 



df 
dq 



T) 



n=0, q=q rBa .i 



Or) 







1=0, q=q rc *i 



which gives as the determining factor 



-q 3 v 3 



V ,|(4sm 2 §-i;V) 



(36) 



(37) 



q=q rC ai 



This condition is exactly what appears in Slepyan's work. 



IV. NUMERICAL CALCULATIONS 

In this section we derive an integral form for the basic result Eq. ([32]). We then present a numerical procedure for 
finding the displacement- velocity curve. The basic notion is to utilize Eq. (|6^) to find S + . It is more convenient for 
the numerical work to isolate explicitly the singularity at q = 0, defining K(q) by 

q z + cp z 

where <f> is an arbitrary positive constant. Because of the factor A 2 q 2 in the denominator, K(0) — 1; the factor 
(q 2 + 4> 2 )/(j) 2 does not alter the desired asymptotic behavior at infinity: K — > const as q — > ±00. Applying Eq. ([63]), 
we get for Ira q > 

s+ _ exp _L ! + °° iifiMj = (^i±«\ \ K+)i m 

^j-oo 1-q V ? + W V ^ 
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and for Im q < 



with 



* =e *»ML ~^ dl (41) 

where the upper sign is for im g > and the lower sign for Im q < 0. 

There is one point regarding these expressions that must be clarified. We must make the correct choice of the branch 
cuts for the square roots in h and r appearing in S. It is convenient to express things in terms of S*(q) = h/r = H/R, 
the simple k = 1 and rjk — f], form of S, so that from Eqn. ( |l5| ) 

5 = S*(l + l V vq) 

S*(l + i7p;g) + (l-S*)fc(l+i77 fc i;g) ' V ' 

and from Eqn. 

1 — <?* 

« = TT^ (43) 

We must choose S* such that |£| < 1 as noted below Eqn. ([)]). It is easy to see from the above expression for £ that 
it is sufficient to choose S* to lie in the first and fourth quadrants to guarantee this. For S*, we can then generate S 
directly via Eqn. (^2), and from this K. 

As K(0) = 1, there is no pole in the integral Q) at 7 = when q = 0. This means that K + (0) = (if "(O)) -1 and 
the asymptotic expressions for 5 ± are 



S + \ q ^ = ] l^r^(q + ^: S~\^ Q = J -^{q - i0)* . (44) 



Comparing this to Eqn. (p3|), we have 



so that 



k = V-i^ + (0) , (45) 



_A l + <fn, k v K+m 



Explicitly writing out the integrals, we get 



A l + <hkv 1 f +co lnif( 7 ) 

ex P7— / I ^ 7 (47) 



A G y ^40 47ri J.^ 7(1 + 477^^7) 

The extra degree of convergence of the integrand at infinity is very desirable, as it helps control the fact that lnif (£) 
is a rapidly oscillating function. Also, the fact that the answer must be independent of <fi provides a check on the 
numerical routines. A change of variables 7 — > —7 allows the integral to be rewritten as 




(48) 

To actually compute this integral, we divided the region of integration into two parts from to 1 and from 1 to +00. 
The second integral was further transformed so as to have in the numerator the factor In A 2 (f> 2 K((,) which goes to 
zero at infinity; the subtracted integral with the numerator lnl/(A<fi) 2 was performed analytically. After all these 
manipulations, the integrals were successfully computed by standard mathematical library subroutines. 

Results from our numerical calculations for k = 1 and rjk = fj = T], are presented in Figs, [l], pi which reproduce, 



of course, those of (Kessler, 2000). At small 77, there are large oscillations in the v vs A curve, as found initially by 
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A/A, 



FIG. 1 Dependence of speed v on dimensionless driving ^— for k = 1 and r\u — r\. The curves are for r\ — 0.01, 0.1, 0.7, 1.3, 
1.9 G 



0.2 



0.1 



n=0.01 
71=0.04 
11-1.7 



!/ ; 
i 

if i 



0.0 L 



1.40 



1.50 



A/A, 



FIG. 2 Change from rapid oscillations for small rj to smooth behavior at large rj. Here k = 1 and r]k = V- 



Marder and Gross for lattices with small transverse sizes. The curves all hit th e v = line at the point which marks 
the end of the band of lattice-trapped static cracks (Kessler and Lcvinc, 1998). For larger damping, the oscillations 
disappear and the moving crack branch bifurcates smoothly (in the backward direction) from this point, eventually 
turning around and becoming stable. In Fig. |^, we contrast the standard result to that for smaller fc, keepin g rjk = r\. 
We see that the v = limit has moved closer to unity, in accord with the finding in ( Kessler and Lcvinc, 1998 ) that the 
window of arrested cr acks narrows as k decreas es (note for purposes of comparison that our k is a factor of 2 smaller 
than that defined in ( Kessler and Levine, 199£ ). The velocity at the minimum A/Ag however is hardly affected, so 
the minimum stable velocity is essentially unchanged. Past this point, the velocity is higher at smaller k, even when 
the change in Aq is factored in. This is further evidence of how the microscopic details influence the velocity versus 
driving displacement curve. 



V. CONSISTENCY OF THE SOLUTION 

As mentioned in the introduction, cracks become unstable above a critical velocity. In the framework of the model 
being considered here, one might assume that this velocity is associated with the solution as found by the Wiener-Hopf 
method becoming inconsistent. As we shall see, there are two kinds of inconsistencies. The first sets in for small v, 
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FIG. 3 Dependence of speed v on dimensionless driving A/Aa for k = 0.01, 0.1, 0.2, 0.4, 1.0. Here r\ — rj k — 1.1 



for small enough r\. Here the vertical bond between y = and y = 1 first achieves critical extension prematurely, at 
some positive r, contradicting the assumption that the bond breaks at r = 0. In the second kind of inconsistency, a 
horizontal bond achieves critical extension. If the model is defined so that only the vertical bonds between y — and 
y = 1 are breakable, then this does not present a problem. If one the other hand, all the springs are chosen identical, 
so that the horizontal springs also break at the same critical extension as the vertical ones, then the solution is indeed 
inconsistent at this point. This is what occurs at large enough velocity, for all r\. 

To proceed, we must calculate the bond lengths, back transforming our Fourier space solutions so as to obtain the 



physical displacements. We find the time dependence of the function V(r). From Eqs. (EOal), (20b) and (felh, we have 



v+ = m s+ w,v(o) , 

(-iq+0)(l + iqt)kv) S+(q)\ a= ^_ l+iq^v' 

v - = V(0) 1 , VkvV(0) m 

S + \ Q= _^S'(q) {iq + 0)(l + iqrikv) l + iqrj k v ' 



Performing the inverse Fourier transform, we get 



A(j)e- %qT (K+)i dq 



V W = V(0) / „. 2. „ 7—. r di ^ (51) 



^{iq-Q){iq -</>)(! + iqr lk v) S+ \ g =^ 



2tt 



for r > and 



«/_«, y/ A<p(iq + 0){1 + iqrjkV) \iq + J S+\ ^_ 2?r 

for r < 0. Here we have used the expressions for derived in the last section. 

Clearly, these integrals must be done numerically. To proceed, we transform the pieces containing the factors K 

1/2 

by dividing both the numerator and denominator by (K + (0)) . In the denominator, we get Ag/A; we have already 
discussed how this can be computed. We now have integrands containing factors of the form 

(^ ± (9)) ±l 1 f + °° q\nK(0 

- ex P7- / 7Ft ~s # ■ ( 53 ) 



The ± refers on the right hand side to the sign of the imaginary part of q. As we will see later on, we need to calculate 
this function only for positive q. Then, the integral in the exponent can be written as 

i r- * >» *«> ,« + • p J^*'*L« . (54) 



2t« Jo i 2 ~ (1 ± «0) 2 s 2tt J £(£ 2 - (1 ± i0) 2 ) 
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We break up the region of integrations into three parts, (0, 1/2), (1/2,3/2) and (3/2, +oo). In the first interval, we 
numerically calculate the integral directly as written. For the third interval, we proceed as discussed in the last section 
for a similar semi-infinite integral; this leaves us with having to integrate numerically a function which behaves as 
l/£ near infinity. Finally, for the integral near 1 we add an integral with K(q£) replaced by K(q); this added integral 
is done analytically and the resultant subtracted integral with the integrand now containing \nK{q^)/K{q) is done 
numerically. 

i i i 

Applying the described procedure, we can evaluate the functions (K ((?)) 2 / (K + (0)) 2 for any real positive q. 
We must now perform the final integration over the variable q. First, we note that the requirement of having a real 
displacement necessitates V ± (— q) = V ± (q)*. This allows us to change the region of integration to (0, +oo) in (|5l|); 



we will return shortly to 52. The trickiest part of this calculation concerns the behavior near q = 0. For V + , there 
is a 1/q 1 / 2 behavior. Here the leading order term can be integrated analytically over the interval (0, 1) and then the 
overall function with this term subtracted can be used for a numerical integration. This works in a straightforward 
manner. 

We need to be more careful with the integral (|52]), because here the integrand behaves near zero as 1/q 3 ^ 2 . To 
proceed, we first divide the interval of integration into three parts: (— oo, — 1), (—1,1), and (l,+oo). The integrals 
over (— oo, —1) and (1, +oo) can again be combined to yield an integral twice the real part of the integrand over the 
latter range. For the range spanning zero, we first subtract from the integrand the leading term 

1 



AS+\^ {iq + 0) 3 / 2 

After this subtraction, the integrand becomes of order 1/q 1 / 2 . Now, the one-sided integral converges, and thus we 
can again transform the range to (0,1). This subtracted integral is evaluated numerically from a very small value 
of q (typically q = 10~ 3 ) to 1. In the remaining part of the range, i.e. from to 10~ 3 , we evaluate the integral 
analytically after approximating it by its leading small q behavior, a/q 1 / 2 ] the value of the constant a is determined 
by fitting the the behavior of the function near the point 10 -3 . Finally we need to perform analytically the integral 
of the subtracted piece l/(\/A S + \ _L_)/{iq + 0) 3 / 2 . To accomplish this, we deform the contour of integration from 

1 7J k V 

(—1, 1) on the real axis to the curve from —1 to 1 going over the lower half of the unit circle. Note that the branch-cut 
for this function must be taken as before over negative real axis; this choice guarantees the fact that f(—q) = /*(?)• 
Note that once all our transformations are complete, we only need values for the integrand for positive (real) values 
of q; this was used earlier in the technique for calculating . A good check of this complex numerical technique of 
computing the Fourier transformation is applying it to the Fourier transform V + (q) for r < and V~(c/) for r > 0, 
as in this case result must be equal to zero. 

We first examine what happens in the standard case, k = 1 and rj^ = i) = r/. Figure ^ shows the behavior of V(r) 
for 7/ = .2 and various values of the velocity. As discussed above, we see the two different kinds of inconsistencies 
exhibited by these solutions. For example, the solution for v = 0.2 is inconsistent as V(r)/ V (0) > 1 for r > 0. This 
type of inconsistency, seen at small velocities, corresponds to the region of backward behavior on v — A graph. As 
already claimed by Marder and Gross, this region is therefore unphysical. 

More interesting is the inconsistency which sets in at large velocity, and has been argued to be related to the 
experimentally observed microbranching. What happens here is that the elongation of the horizontal bond, given by 
V/iot-M = (V(t — 1) — V(r))/2 increases as a function of increasing driving. At some critical speed it reaches the value 
V(0) and the bond breaks. Figure || shows the dependence of the critical speed on 77. We note in passing that these 
graphs agree with what can be gotten by extrapolating a s imilar graph for the finite N y lattice, a graph that can be 



obtained using the methods of ( |Kessler and Levine, 1998 ). The actual location in r at which the horizontal break 
occurs at the critical velocity is shown in Fig. |6[ 

Several details of our findings are worthy of note. First, there is always a critical speed for any dissipation, but this 
speed gets very close to the maximal crack speed for large 77. The value of the critical speed for the dissipation-less 
limit is v cr ~ .725; this is above what has been seen in some experiments but recall that we are doing mode III on a 
square lattice, a far cry from mode I in an amorphous system. The is a curious break in the curve at 77 ~ .665, which 
corresponds in Fig. || to the point where the break location varies most strongly with 7/. 



It should be noted that this inconsistency does not appear related in any obvious way to the Yoffe (1951 ) criterion 



First, the Yoffe criterion, concerning the direction of maximal stress, is calcula ted for the continuum elastic field, 



and is thus completely independent of the viscosity 77 ( Kcsslcr and Levine, 1998 ). The onset of the inconsistency is 



strongly 77 dependent however. Also, the Yoffe criterion is stated for Mode I cracks, and its applicability to Mode III 
is subject to question. 

We turn now to the case of k < 1, with weakened bonds along the crack path, again keeping 77^ = 7/. Figures 
H H show our results. As is reasonable, the inconsistencies persist, but are pushed to higher speeds as k decreases. 
It is interesting to speculate that since for finite width systems, velocities greater than the wave speed are possible, 
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FIG. 4 Elongation of vertical bonds along the path of the crack for r] = 0.2. Here k = 1 and rjk = r\ 
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FIG. 5 Dependence of the critical speed v cr on rj for k = 1, 0.75, 0.5. Here r/k = 77 



sufficiently small k might enable stable steady-state fracture at supersonic velocities. Also, intersonic wave speeds are 
possible in the case of mode II fracture, and therefore weakened bonds on the path of the crack would in this case as 
well allow steady-state propagation at higher speeds (Rosakis, et al., 1999). 



VI. DISCUSSION 



We have shown how to extend the Slepyan approach to cracks in ideally brittle infinite lattice systems to the case 
of having a Kelvin viscosity for each lattice spring. In addition, we have examined the onset of the additional, off-axis 
cracking at high velocity. We showed how dissipation delays but does not eliminate the microbranching instability, 
an effect impossible to see within the context of any continuum elastic treatment. 

In the following work, we will present the results of similar calculations for mode I cracks on a triangular lattice. 
That system is much closer to those that have been studied experimentally and one can therefore hope for more 
direct lessons to emerge. Also, we should mention complementary studies on nonlinear lattice models which relax the 
assumption of ideally brittle springs and thereby allow one to view the onset of inconsistent behavior as the limiting 
case of more standard bifurcations (Kcsslcr and Levinc, 1999, 2000| ). 
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FIG. 6 Dependence of the critical time r cr on r\ for k = 1, 0.75, 0.5. Here rjk = r\. 
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VII. APPENDIX 

We discuss here some of our conventions and notations regarding Fourier transformations. We use Fourier trans- 
formation in the form 

/+oo 
u(T)e lqT dT 
-OO 

i r+oc 

u(r) = —J u F (q)e-^dq (55) 

We define also u + and u~ as 



This gives 



u + (q) = / M(r)e l9T dr (56a) 
Jo 

u-(q) = { u{r)e lqT dT (56b) 

J — OO 

u{t)9{t) - — / u+(q)e-^dq ; (57a) 



— OO 



u{t)6{-t) = -!- / u-{q)e- lqT dq , (57b) 

27T 



which means that u + has poles only in the lower half plane and u only in the upper half plane. 
Now let us define p — —iq. Then 

u{r)e lqT dT P ~°° u(0) / e- pT dr = -«(0) . (58) 
Jo P 
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Thus 



Similarly we can get for p — iq 



So 



lim -iqu + = u(0) . (59) 

-iq^-\-oo 



u(T)e lqT dT P ~ u(0) / e pr dr = -u(0) . (60) 



lim iqu = u(0) . (61) 



The Fourier transform of Q(—t)-^u{t) is given by 

p0 du{r) 
dr 



f 



e lqT dr = u(0) - iqu' (62) 



Finally, we consider separating an arbitrary function S(q) into the product of two pieces S + and S , with poles and 
zeroes in lower and upper half-planes respectively. These can be done using the identity, 



for Ira q > 0, lnS"^) 
/or Jm g < 0, — In S~ (q) 
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